{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "81d831b3b0e92996",
   "metadata": {},
   "source": [
    "### 基于遗传编程自动设计优化算法\n",
    "众所周知，演化计算中一个重要的研究课题就是设计新的优化算法。这个过程通常是由人类专家完成的，但是，我们是否可以让计算机自动设计优化算法呢？这个问题的答案是肯定的。本文将介绍如何基于遗传编程自动设计优化算法。\n",
    "\n",
    "**根据这样一个自动算法设计的工具，我们在得到一个算法公式之后，只要再观察一下自然界中是否有对应的生物行为，就可以得到一个新的智能优化算法。**\n",
    "\n",
    "比如，本文将尝试使用遗传编程自动设计出北极狐算法！\n",
    "\n",
    "![北极狐算法](img/Fox2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e3427b91e831fc9",
   "metadata": {},
   "source": [
    "### 优化函数\n",
    "比如，我们希望自动设计出的算法可以再球型函数上表现良好。球型函数是一个单目标优化领域中的经典测试函数，其公式如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "initial_id",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-02-07T23:56:31.688305600Z",
     "start_time": "2024-02-07T23:56:31.666788Z"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import operator\n",
    "import random\n",
    "\n",
    "from deap import base, creator, tools, gp, algorithms\n",
    "import numpy as np\n",
    "\n",
    "np.random.seed(0)\n",
    "random.seed(0)\n",
    "\n",
    "\n",
    "def sphere(x, c=[1, 1, 1]):\n",
    "    \"\"\"\n",
    "    Shifted Sphere function.\n",
    "\n",
    "    Parameters:\n",
    "    - x: Input vector.\n",
    "    - c: Shift vector indicating the new optimal location.\n",
    "\n",
    "    Returns:\n",
    "    - The value of the shifted Sphere function at x.\n",
    "    \"\"\"\n",
    "    return sum((xi - ci) ** 2 for xi, ci in zip(x, c))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d87e994c3144076d",
   "metadata": {},
   "source": [
    "### 经典优化算法\n",
    "在文献中，差分演化可以用来求解这个球型函数优化问题。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "feb772104d562277",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-02-07T23:56:31.817414Z",
     "start_time": "2024-02-07T23:56:31.695306200Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "传统DE算法得到的优化结果 4.506377260849465e-05\n"
     ]
    }
   ],
   "source": [
    "# DE\n",
    "dim = 3\n",
    "bounds = np.array([[-5, 5]] * dim)\n",
    "\n",
    "\n",
    "# Define a simple DE algorithm to test the crossover\n",
    "def differential_evolution(\n",
    "        crossover_func, bounds, population_size=10, max_generations=50\n",
    "):\n",
    "    population = [\n",
    "        np.random.rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0]) + bounds[:, 0]\n",
    "        for _ in range(population_size)\n",
    "    ]\n",
    "    population = np.array(population)\n",
    "    best = min(population, key=lambda ind: sphere(ind))\n",
    "    for gen in range(max_generations):\n",
    "        for i, x in enumerate(population):\n",
    "            a, b, c = population[np.random.choice(len(population), 3, replace=False)]\n",
    "            mutant = np.clip(crossover_func(a, b, c, np.random.randn(dim)), bounds[:, 0], bounds[:, 1])\n",
    "            if sphere(mutant) < sphere(x):\n",
    "                population[i] = mutant\n",
    "                if sphere(mutant) < sphere(best):\n",
    "                    best = mutant\n",
    "    return sphere(best)\n",
    "\n",
    "\n",
    "print(\"传统DE算法得到的优化结果\",\n",
    "      np.mean([differential_evolution(lambda a, b, c, F: a + F * (b - c), bounds) for _ in range(10)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e46b8aec8871cdd9",
   "metadata": {},
   "source": [
    "可以看到，传统DE算法得到的优化结果是不错的。但是，我们是否可以自动设计出一个更好的算法呢？"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "712f8d2a7147ff03",
   "metadata": {},
   "source": [
    "### 基于遗传编程的自动设计优化算法\n",
    "其实DE的交叉算子本质上就是输入三个向量和一个随机向量，然后输出一个向量的函数。因此，我们可以使用遗传编程来自动设计这个交叉算子。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "3b598a4e994266e8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-02-07T23:56:46.285724800Z",
     "start_time": "2024-02-07T23:56:31.818414300Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gen\tnevals\tavg   \tmin      \tmax    \n",
      "0  \t50    \t2.6796\t0.0112234\t15.2248\n",
      "1  \t50    \t2.41407\t0.00253387\t17.9657\n",
      "2  \t45    \t1.41727\t0.0205569 \t18.5921\n",
      "3  \t47    \t0.99445\t0.00658522\t14.4601\n",
      "4  \t47    \t0.929668\t0.005623  \t13.84  \n",
      "5  \t48    \t1.61888 \t0.00913134\t13.9251\n",
      "6  \t50    \t1.18172 \t0.000383948\t14.9727\n",
      "7  \t48    \t0.624159\t0.000705421\t12.3018\n",
      "8  \t50    \t0.765903\t0.00214913 \t8.71667\n",
      "9  \t43    \t0.3652  \t0.0110385  \t3.56652\n",
      "10 \t47    \t1.39889 \t0.00685267 \t22.123 \n",
      "11 \t43    \t1.27877 \t0.00685267 \t20.31  \n",
      "12 \t48    \t1.82377 \t0.0027862  \t11.4693\n",
      "13 \t49    \t0.736725\t0.0108848  \t12.7022\n",
      "14 \t50    \t1.39344 \t0.0102804  \t12.8329\n",
      "15 \t47    \t0.847688\t0.00398283 \t11.3424\n",
      "16 \t44    \t0.9867  \t0.0067096  \t15.8511\n",
      "17 \t48    \t0.971622\t0.0180985  \t9.05041\n",
      "18 \t42    \t0.843393\t0.00948021 \t11.9563\n",
      "19 \t47    \t0.849741\t0.00759852 \t10.9686\n",
      "20 \t47    \t0.999861\t0.00425035 \t14.4111\n",
      "21 \t42    \t1.18842 \t0.00665311 \t13.5106\n",
      "22 \t46    \t1.41895 \t0.00320289 \t15.9007\n",
      "23 \t47    \t1.19332 \t0.00406941 \t9.579  \n",
      "24 \t48    \t0.923953\t0.00313277 \t11.4326\n",
      "25 \t45    \t0.599486\t0.00469191 \t8.87691\n",
      "26 \t43    \t1.06541 \t3.39457e-29\t15.4452\n",
      "27 \t44    \t1.38335 \t0.00224764 \t13.3298\n",
      "28 \t48    \t1.45239 \t0.017065   \t9.51407\n",
      "29 \t48    \t1.08886 \t0.00518668 \t12.8216\n",
      "30 \t48    \t0.55234 \t0.00209358 \t6.49766\n",
      "Best Crossover Operator:\n",
      "add(ARG0, subtract(multiply(ARG0, ARG3), ARG3))\n",
      "Fitness: (3.3945670827791664e-29,)\n"
     ]
    }
   ],
   "source": [
    "# GP 算子\n",
    "pset = gp.PrimitiveSetTyped(\"MAIN\", [np.ndarray, np.ndarray, np.ndarray, np.ndarray], np.ndarray)\n",
    "pset.addPrimitive(np.add, [np.ndarray, np.ndarray], np.ndarray)\n",
    "pset.addPrimitive(np.subtract, [np.ndarray, np.ndarray], np.ndarray)\n",
    "pset.addPrimitive(np.multiply, [np.ndarray, np.ndarray], np.ndarray)\n",
    "pset.addEphemeralConstant(\"rand100\", lambda: np.random.randn(dim), np.ndarray)\n",
    "\n",
    "pset.context[\"array\"] = np.array\n",
    "\n",
    "creator.create(\"FitnessMin\", base.Fitness, weights=(-1.0,))\n",
    "creator.create(\"Individual\", gp.PrimitiveTree, fitness=creator.FitnessMin)\n",
    "\n",
    "toolbox = base.Toolbox()\n",
    "toolbox.register(\"expr\", gp.genHalfAndHalf, pset=pset, min_=1, max_=2)\n",
    "toolbox.register(\"individual\", tools.initIterate, creator.Individual, toolbox.expr)\n",
    "toolbox.register(\"population\", tools.initRepeat, list, toolbox.individual)\n",
    "toolbox.register(\"compile\", gp.compile, pset=pset)\n",
    "\n",
    "\n",
    "# Evaluate function for GP individuals\n",
    "def evalCrossover(individual):\n",
    "    # Convert the individual into a function\n",
    "    func = toolbox.compile(expr=individual)\n",
    "    return (differential_evolution(func, bounds),)\n",
    "\n",
    "\n",
    "toolbox.register(\"evaluate\", evalCrossover)\n",
    "toolbox.register(\"select\", tools.selTournament, tournsize=3)\n",
    "toolbox.register(\"mate\", gp.cxOnePoint)\n",
    "toolbox.register(\"mutate\", gp.mutUniform, expr=toolbox.expr, pset=pset)\n",
    "\n",
    "# Evolve crossover strategies\n",
    "population = toolbox.population(n=50)\n",
    "hof = tools.HallOfFame(1)\n",
    "stats = tools.Statistics(lambda ind: ind.fitness.values)\n",
    "stats.register(\"avg\", np.mean)\n",
    "stats.register(\"min\", np.min)\n",
    "stats.register(\"max\", np.max)\n",
    "\n",
    "algorithms.eaSimple(population, toolbox, 0.9, 0.1, 30, stats, halloffame=hof)\n",
    "\n",
    "# Best crossover operator\n",
    "best_crossover = hof[0]\n",
    "print(f\"Best Crossover Operator:\\n{best_crossover}\")\n",
    "print(f\"Fitness: {best_crossover.fitness.values}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf2c8fb3e5d148c6",
   "metadata": {},
   "source": [
    "### 分析新算法\n",
    "现在，我们得到了一个新的交叉算子。我们可以看一下这个交叉算子的公式。\n",
    "$X_{new}=X+(F*X-F)$, F是一个随机变量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "id": "71c1e9de586767b0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-02-07T23:58:03.859051200Z",
     "start_time": "2024-02-07T23:58:03.730618Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "新优化算法得到的优化结果 1.0213225557390857e-19\n"
     ]
    }
   ],
   "source": [
    "add = np.add\n",
    "subtract = np.subtract\n",
    "multiply = np.multiply\n",
    "square = np.square\n",
    "array = np.array\n",
    "\n",
    "crossover_operator = lambda ARG0, ARG1, ARG2, ARG3: add(ARG0, subtract(multiply(ARG0, ARG3), ARG3))\n",
    "print(\"新优化算法得到的优化结果\", np.mean([differential_evolution(crossover_operator, bounds) for _ in range(10)]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c39ad9e7553bc87",
   "metadata": {},
   "source": [
    "从结果可以看到，新的优化算法得到的优化结果优于传统DE算法。这证明GP发现了一个更好的新算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37941230eb02cbab",
   "metadata": {},
   "source": [
    "### 北极狐算法\n",
    "现在，这个算法我们可以命名为北极狐算法。北极狐的毛色会根据季节变化。在这个公式中，X会根据随机变量F的变化而变化。这个公式的形式和北极狐的毛色变化有些相似。因此，我们可以将这个算法命名为北极狐算法。\n",
    "![北极狐算法](img/Fox.png)\n",
    "\n",
    "该算法的交叉算子为$X_{new}=X+(F*X-F)$。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
